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ABSTRACT 



The orbital separation of compact binary stars will shrink with time due to the 
emission of gravitational radiation. This inspiralling phase of a binary system's 
evolution generally will be very long compared to the system's orbital period, but the 
. final coalescence may be dynamical and driven to a large degree by hydrodynamic 

qq , effects, particularly if there is a critical separation at which the system becomes 

dynamically unstable toward merger. Indeed, if weakly relativistic systems (such 
| as white dwarf-white dwarf binaries) encounter a point of dynamical instability at 

some critically close separation, coalescence may be entirely a classical, hydrodynamic 
process. Therefore, a proper investigation of this stage of binary evolution must include 
three-dimensional hydrodynamic simulations. 
^ ■ We have constructed equilibrium sequences of synchronously rotating, equal-mass 

r*"** , binaries in circular orbit with a single parameter - the binary separation - varying 

along each sequence. Sequences have been constructed with various polytropic as 
q-i well as realistic white dwarf and neutron star equations of state. Using a Newtonian, 

5^ | finite-difference hydrodynamics code, we have examined the dynamical stability of 

individual models along these equilibrium sequences. Our simulations indicate that 
no points of instability exist on the sequences we analyzed that had relatively soft 
^ ■ equations of state (polytropic sequences with polytropic index n = 1.0 and 1.5 and 

two white dwarf sequences). However, we did identify dynamically unstable binary 
models on sequences with stiffer equations of state (n = 0.5 polytropic sequence and 
two neutron star sequences). We thus infer that binary systems with soft equations of 
state are not driven to merger by a dynamical instability. For the n = 0.5 polytropic 
sequence, the separation at which a dynamical instability sets in appears to be 
associated with the minimum energy and angular momentum configuration along the 
sequence. Our simulations suggest but do not conclusively demonstrate that, in the 
absence of relativistic effects, this same association may also hold for binary neutron 
star systems. 



Subject headings: binaries: close — hydrodynamics — instabilities — stars: neutron 
— white dwarfs 
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1. Introduction 

The coalescence of double white dwarf and double neutron star binaries is important to 
examine since this process may produce a number of astrophysically interesting objects and 
events. Double white dwarf binary mergers have been suggested as precursors to some Type la 
Supernovae (Iben & Tutukov 1984; Iben 1988, 1991; Iben & Webbink 1989; Branch et al. 1995) 
and to long gamma-ray bursts (Katz & Canel 1996). White dwarf-white dwarf mergers may also 
lead to the formation of massive single white dwarfs or neutron stars (Colgate & Petschek 1982; 
Saio & Nomoto 1985; Iben & Tutukov 1986; Kawai, Saio, & Nomoto 1987; Chen & Leonard 1993); 
to the formation of subdwarf stars; or to the formation of hydrogen deficient, highly luminous 
stars (Iben 1990 and references therein; Webbink 1984). Neutron star-neutron star mergers may 
result in bursts of gamma-rays and neutrinos, in the production of r-process elements, and in the 
formation of black holes (Eichler et al. 1989; Meyer 1989; Narayan, Paczhski, & Piran 1992; Rasio 
& Shapiro 1992; Davies, Benz, Piran, & Thielemann 1994; Katz & Canel 1996; Lipunov et al. 
1995; Ruffert, Janka, h Schafer 1996; but see the simulations of Shibata, Nakamura, and Oohara 
[1993] and Janka & Ruffert [1996] which cast doubt on the neutron star-neutron star merger 
scenario as a precursor to gamma-ray bursts). 

Merging compact binaries are also expected to be relatively strong sources of gravitational 
radiation. The gravitational radiation emitted during the inspiral phase of double neutron star 
binary evolution (i.e., before tidal effects become sizeable) is likely to be detected by terrestrial 
interferometric detectors such as LIGO and VIRGO, which will be sensitive to frequencies in 
the range of 10-10 3 Hz (Abramovici et al. 1992; Cutler et al. 1993; Thorne 1995). Proposed 
"dual-recycled" interferometers and spherical "TIGA" type resonant detectors will be more 
sensitive than LIGO to the higher frequency radiation, ~ 10 3 Hz, emitted during the brief final 
merger stage of the coalescence (Meers 1988; Strain & Meers 1991; Cutler et al. 1993; Merkowitz 
& Johnson 1994; Thorne 1995; however, see Wilson & Mathews 1995 who indicate that the 
gravitational wave radiation emitted during this phase may have a lower frequency than previously 
expected). The gravitational wave radiation emitted during the merger phase in double white 
dwarf binary evolution is unlikely to be detected in the near future because the expected frequency 
of the radiation falls in or just beyond the upper end of the frequency range (10 _4 -10 _1 Hz) of 
proposed space-based laser interferometric detectors (Faller et al. 1989; Hough et al. 1995) and 
the duration of the phase will be too short to produce a significant signal in this range of the 
detectors' sensitivity. 

Because the final stages of binary coalescence are driven in part by sizeable nonlinear tidal 
effects, numerical hydrodynamic techniques must be used to properly follow the evolution of 
merging binaries. The first step in performing such a hydrodynamic simulation is the construction 
of an appropriate initial model. The coalescence of the chosen binary system must proceed on 
a dynamical timescale (on the order of a few initial orbital periods), in order for an explicit 
hydrodynamics code to be able to carry out the simulation in a reasonable amount of computational 
time. Hence, the components of the initial binary model must either be at a separation where they 
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are dynamically unstable to coalescence or they must be forcibly brought to coalescence from a 
wide separation (e.g., by draining orbital angular momentum away from the system in a way that 
mimics the effects of the gravitational wave radiation reaction). Using the former methodology, 
the work presented herein focuses on the identification of dynamically unstable binary systems. 

The initial separation at which a particular binary model becomes dynamically unstable to 
merger, if one exists, can be found via a stability analysis of a set of binary models constructed 
in hydrostatic equilibrium, along a constant mass sequence of decreasing orbital separation. This 
sequence serves as an approximate representation of the evolution of the binary as its components 
are brought closer together by the effects of gravitational radiation. Such analyses have recently 
been done by Lai, Rasio, & Shapiro (hereafter LRS) and by Rasio & Shapiro (hereafter RS) for 
binaries with polytropic equations of state (EOS). In a polytropic equation of state, the pressure 
P is expressed in terms of the density p as P = Kp 1+1 / n , where K is the polytropic constant 
and n is the polytropic index (see §2). The analytical work of LRS utilized an approximate 
energy variational method and studied detached binaries with components having various mass 
ratios, spins, and polytropic indices (LRS 1993a, b, 1994a, b). The numerical work of RS utilized 
the smoothed particle hydrodynamics technique to study detached and contact binaries with 
components having various mass ratios, but equal spins and polytropic indices (RS 1992, 1994, 
1995; for earlier work see Gingold & Monaghan 1979; Hachisu & Eriguchi 1984a, b; Hachisu 
1986b). 

We performed stability analyses of equilibrium sequences of double white dwarf binaries 
constructed with the zero-temperature white dwarf equation of state (Chandrasekhar 1967), 
double neutron star binaries constructed with realistic neutron star equations of state (adapted 
from Cook, Shapiro, & Teukolsky 1994), and, for the sake of comparison with the work of LRS and 
RS, polytropic binaries with n = 0.5, 1.0, and 1.5 equations of state. The examined equilibrium 
sequences were constructed with the Self-Consistent Field technique of Hachisu (1986a, b), 
which produces models of rotating, self-gravitating fluid systems in hydrostatic equilibrium. For 
simplicity, all binary models along these sequences were constructed as synchronously rotating 
systems having equal-mass (q = 1.0) components. The relative stability of individual binary 
systems along selected sequences was examined using a three-dimensional (3D), finite-difference 
hydrodynamics code. Both the construction of our equilibrium binary sequences and our stability 
tests along these sequences have been done using purely Newtonian gravity and Newtonian 
dynamics. 

Our numerical techniques are briefly described in Section 2. Constructed equilibrium 
sequences are presented in Section 3 and our dynamical tests of the stability of individual models 
along selected sequences are presented in Section 4. Finally, the implications of these results are 
discussed in Section 5. 
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2. Numerical Techniques 

Our simulations of close binary systems involve the solution of the following set of equations 
which govern the structure and evolution of a nonrelativistic fluid in cylindrical coordinates: 

the continuity equation, 

|? + V.(H = 0; (1) 
the three components of the equation of motion, 

9S „ dP d<$> A 2 

^ + V. W = ---„- + — 3 , (2) 
dT „ . dP d<S> 

_ + V .(T l0 = (3) 
_ + V .(A„) = Si-Pgr; (4) 



Poisson's equation 



V 2 $ = AttG P ; (5) 



and the equation of state (see below). In the above equations, v is the velocity; S = pu, T = pw, 
and A = pRv<p are the radial, vertical, and angular momentum densities, respectively (where u, w, 
and are the radial, vertical, and azimuthal components of the velocity, respectively); R, 4>, and 
z are the cylindrical coordinates; and $ is the gravitational potential. 

We have used three types of barotropic equations of state in this work. The first, and simplest, 
type is a polytropic equation of state for which 

P = Kp l+l ' n , (6) 

where K is the polytropic constant and n, the polytropic index, determines the degree of 
compressibility of the fluid (the higher the value of n, the more compressible/softer the fluid). 

The second type of equation of state used is the zero-temperature white dwarf (WD) 
equation of state (Chandrasekhar 1967), which represents the pressure distribution of a completely 
degenerate electron gas: 



P = a 



i/o _ 

x (2x 2 - 3) (x 2 + l) + 3 In (x + Vl + x 2 ) 



X EE (p/b f 3 , (7) 

where ao = 6.00 x 10 22 dynes cm" 2 , 60 = l-95(// e /2) x 10 6 gcm -3 , and p e is the mean molecular 
weight per electron. We have used p e = 2 in all of our computations. The heaviest nonrotating 
single object that can be constructed with this equation of state has a mass of 1.44M© (this is the 
Chandrasekhar mass Men)- 
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The third type of equation of state used here is a realistic neutron star (NS) equation of 
state. We have chosen three such equations of state (from among the fourteen realistic NS 
equations of state listed in Cook, Shapiro, & Teukolsky [1994, hereafter CST]), each with a 
different degree of compressibility (one soft, one medium, and one hard). Specifically, the chosen 
soft equation of state is CST's equation of state F; the medium one is CST's equation of state 
FPS; and the hard one is CST's equation of state L (see references within CST for the original 
sources of these equations of state). We obtained these equations of state in tabular form from 
Cook (1995). The tables each provide ~ 500 values of the pressure P for values of p ranging 
over 15 orders of magnitude, from ~ 8 g cm -3 to 10 16 g cm -3 (note that it is actually the number 
density N = p/m neutr0 n-, where m neu t r on = 1-67 x 10~ 24 g, that is tabulated). Because we wanted 
to perform parallel finite-difference hydrodynamics (FDH) simulations of systems with these 
equations of state and did not possess an interpolation algorithm designed for efficient use on a 
parallel machine, polynomial fits to the tabular data were necessary. Some numerical manipulation 
of the data was also needed because of the particulars of the technique used in the initial model 
construction. (See New [1996] for details.) 

If the only motion of a fluid system is rotation about an axis with an angular velocity f2, 
which is constant in time and a function only of the distance from the rotation axis, the structure 
of the system is described by the following single expression: 

-VP + V$ + W(i?) = 0, (8) 

where the z-axis has been chosen as the axis of rotation and the centrifugal potential 
ty(R) = — J Q 2 (R)RdR. Such a fluid is said to be in hydrostatic equilibrium because the forces 
due to its pressure and to its gravitational and centrifugal potentials are in balance. All of the 
initial equilibrium binary systems studied in this work have been constructed in hydrostatic 
equilibrium according to this prescription, along with the additional constraint that angular 
velocity is a spatial constant f^o (i- e -) n °t a function of R). In this case of uniform rotation, 

#(#) = -n 2 R 2 /2. 



2.1. Self Consistent Field Code 

The method we have used to construct the equilibrium models is Hachisu's grid-based 3D, 
Self-Consistent Field (HSCF) technique (Hachisu 1986a, b). This iterative technique produces 
rotating, self-gravitating fluid systems in hydrostatic equilibrium. Our version of the HSCF 3D 
code computes the gravitational potential via a direct numerical solution of Poisson's equation 
(5). Details of the method used can be found in Tohline (1978). 

An estimate of the quality of the converged equilibrium configuration is obtained from a 
determination of how well the energy is balanced in the system. This balance is measured by the 
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virial error VE: 

VEs \w±™+m£tl, (9) 

where T is the kinetic energy; W is the gravitational potential energy; and V is the volume of 
the model. The virial errors in our equilibrium models constructed with polytropic and WD EOS 
were typically ~ 10 _3 -10 -4 ; those in models constructed with the realistic NS EOS were typically 
~ 1(T 2 . 

The forms of the WD and realistic NS equations of state are such that when they are used 
in the HSCF code, the density maxima p ma x of the models to be constructed must be given to 
the code as input. Thus because we were interested in constant mass sequences for our stability 
analyses of close binaries, we had to perform an iteration in the choice of p m ax until we arrived at 
a configuration with the desired Mt, in the case of models with the WD and realistic NS equations 
of state. However, in the polytropic case, converged models can actually be obtained without an a 
priori choice of p ma x and then later scaled as desired. 

Our 3D equilibrium configurations are assumed to be symmetric about the z = (equatorial) 
plane; this symmetry will be referred to as equatorial symmetry. Because our version of the 3D 
HSCF code constructs binaries with only equal-mass components, a periodic symmetry over the 
azimuthal range < <f> < n is also assumed. This means that a quantity U specified at an angle <p 
is equivalent to that same quantity specified at all angles ft for which ft = (cf> + mir) and m is an 
integer: 

U{(t> + mn) = U{ft). (10) 
This symmetry will be referred to as 7r-symmetry. 



2.2. Finite Difference Hydrodynamics Code 

A finite-difference hydrodynamics (FDH) code was used to solve, on a discrete numerical 
grid, equations (1-5) and an equation of state (see above), which govern the temporal evolution of 
a fluid. FDH codes differ from smoothed particle hydrodynamics (SPH) codes in that they follow 
the evolution of the fluid as it flows through a fixed set of grid cells, instead of treating the fluid 
as a set of particles and following the evolution of each particle. 

The 3D FDH code used in this present study is a Fortran 90 version of the Fortran 77 
code described by Woodward (1992) and Woodward, Tohline, & Hachisu (1994). It was written, 
principally by Woodward, to take advantage of the parallel architecture of the MasPar computers 
on which it is run. The accuracy of the code is second order in both time and space. The 
numerical techniques employed are discussed in detail in Woodward (1992). The solution to 
Poisson's equation (5) is obtained through the ADI method (Cohl et al. 1997). 

As in the HSCF code, the grid cells are uniformly spaced in each of the three directions 
and equatorial and 7r-symmetries are assumed. The single precision hydrodynamic simulations 
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presented here were performed on cylindrical grids with resolutions of 64 x 64 x 64. 

In the binary stability analysis simulations presented in §4, from ~ 2, 600 to 16,700 timesteps 
were required to follow each binary through one initial orbital period Pj. On the 8K-node MasPar 
MP-1 at LSU, each timestep took ss 19 CPU seconds (or « 73 //sec per grid zone); hence, these 
simulations required ~ 14 to 88 CPU hours per Pj. This large range is due to the variation in 
the size of the integration timestep that could be taken in the different simulations. The size of 
this timestep is restricted in order to ensure the numerical stability of the computations. The 
simulations we performed varied in length from 1 Pj to 5 Pi. A few of our binary dynamical 
stability test simulations were run on the 4K-node MasPar MP-2 at the Scalable Computing 
Laboratory of the DOE Ames Laboratory at Iowa State University. The CPU time per timestep 
required for simulations conducted on the MP-2 is about 3/5 that of the time required to run on 
the MP-1 at LSU. 

Our FDH code typically follows the fluid evolution in the inertial reference frame. However, 
we chose to incorporate the option of running the code in a frame of reference which rotates with 
the initial angular velocity of the fluid, Slo- This choice was motivated by a desire to minimize 
numerical effects which might artificially influence the stability of the binary systems studied. 
The particular effect we sought to minimize was dissipation due to numerical viscosity, which 
arises from the coarseness of the finite differencing. The hope was that diminishing the motion 
of the fluid through the grid by running in the rotating reference frame would also diminish the 
dissipative effects of numerical viscosity on the fluid (see §4 for further details). (We also tried 
updating the angular velocity of the rotating frame of reference once during some of the simulations 
presented here, in order to further minimize the dissipation due to numerical viscosity). 

The rotating reference frame adds two terms to the radial equation of motion (2) and one to 
the azimuthal equation (4) (cf., Norman & Wilson 1978): 

dS „ dP d$ A 2 n2n 2fl A , 

w + v.(*o = + — 3 + + (ii) 

f)A BP r)<b 

w + = -_-„_- 2tw (12) 

The p^qR term that has been added to the radial equation of motion results from the centrifugal 
force; the other two added terms result from the Coriolis force. Note that the centrifugal term in 
equation (11) can be rewritten as A' 2 /(pi? 3 ), where A' = pVL^R 2 . We use this form in the actual 
computation of the centrifugal term, with A' centered at the same place in each grid cell as is A, 
in order to be numerically consistent with the computation of the curvature term A 2 /(pR?) in the 
radial equation of motion. 

A discussion of the boundary conditions, vacuum treatment, and rotation axis treatment 
implemented in our hydrodynamics code is presented in New (1996). 
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3. Equilibrium Sequences 

We have constructed hydrostatic equilibrium sequences of synchronized close binaries with 
polytropic as well as realistic white dwarf and neutron star equations of state. The individual 
binary models along each sequence have the same equation of state and constant total mass Mr, 
but decreasing binary separation a. Here a is the distance measured between the pressure [density] 
maxima of the stellar components. Each such sequence represents a quasistatic approximation 
to the evolution of a binary system in which gravitational radiation gradually carries away the 
system's orbital angular momentum. These binary models were constructed, on 128 x 128 x 128 
grids, with the HSCF technique (see §2.1), which creates models of rotating, self-gravitating fluid 
systems in hydrostatic equilibrium. 

It should be noted that the true physical viscosity present in double neutron star binaries is 
not expected to be strong enough to enforce synchronization (Bildsten & Cutler 1992; Kochanek 
1992) and the viscosity of the degenerate material in double white dwarf binaries probably is not 
strong enough to synchronize them either (this can be shown by applying the arguments given 
in Bildsten & Cutler 1992 to white dwarfs and using the values for the viscosity of degenerate 
material given in Durisen 1973). However, if magnetic fields are present, they may bring about 
synchronization. In any case, synchronization is at least a simplifying assumption. Furthermore, 
since neutron stars have relatively strong gravitational fields, Newtonian models and simulations 
of double neutron star binaries need to be viewed with caution. 



3.1. Polytropic Sequences 

The equilibrium sequences that we constructed for binary models with polytropic indices 
n = 0.5, 1.0, and 1.5 are displayed in Figure 1. For each n, the total energy E and total angular 
momentum J are plotted versus the separation a. Note that we do not claim that the results given 
in this figure or in the rest of the manuscript are necessarily accurate to the number of digits in 
which they are reported; the number of digits in which the results are presented allows the display 
of characteristics of the equilibrium sequences and the differentiation between individual models 
on these sequences. As mentioned in §2.1, the equations of state of the polytropic models are such 
that the total mass of the system Mt does not have to be chosen before its construction, but can 
be scaled afterwards as desired. 

There are in fact three parameters, M = l/2My, Rs p hi an d the polytropic constant K, 
which set the scale of a polytropic model. Here Rs P h is the radius of a spherical star of mass 
M and polytropic index n. These parameters are related according to the following equation 
(Chandrasekhar 1967): 



M = Airm, 



'{n + l)K~ 


2(l-n) 


Rsph 


4ttG 







where m n and r n are Lane-Emden constants for a particular value of n (see Table 1 for their values 
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corresponding to n = 0.5, 1.0, and 1.5). Thus, only two of these three parameters are independent; 
if two of them are specified, the other one is automatically determined. The quantities in Figure 1 
are themselves normalized to G, M, and Rs P h- Specifically, E has been divided by (GM 2 Rg^ h ); J 

has been divided by (G 1 / 2 M 3 / 2 i?^); and a has been divided by Rs P h- 

In the discussion that follows, as in Figure 1, all values of the binary separation a have been 
normalized to Rsph- On each sequence, Point C marks the system where the surfaces of the two 
binary components just come into contact. Systems to the right of this point are detached binaries 
and systems to left are contact binaries or "dumbbells." For the sake of illustration, Figure 2 
displays an isodensity surface image of an example detached binary model (a = 3.28, n = 1.0) and 
of an example dumbbell model (a = 2.70, n = 1.0). Points M mark the models along each sequence 
which have the minimum total energy and the minimum total angular momentum. Along all three 
polytropic sequences, the minimum in E occurs at the same separation as the minimum in J. 
(See §3.4 below for a discussion of the significance of these minima.) The model with the smallest 
separation on each sequence, marked with a T, will be referred to as the "terminal" model. No 
synchronously rotating binary models with the equation of state particular to that sequence can 
be constructed in equilibrium with a smaller separation than that of the terminal model because 
the centrifugal force would exceed the gravitational force along the equator of such systems. 

The separation at which this termination of the sequence occurs increases from a = 2.45 for 
n = 1.5 to a = 2.76 for n = 0.5. The separation at which the minima occur also increases from 
a = 2.70 for n = 1.5, to a = 2.89 for n = 1.0, to a = 3.11 for n = 0.5. Note that contact occurs to 
the right of the minima for n = 1.5 and to the left of the minima for n = 1.0 and 0.5. We have 
determined that Points C and M coincide for n = 1.177. If the binary components were spherical, 
their separation at the point of contact would be a = 2. However, a(C) ranges from 2.81 for 
n = 1.5 to 2.94 for n = 0.5. 



3.1.1. Comparison with Previous Work 

In this section, we compare our polytropic equilibrium sequences to those of LRS (1993b) and 
RS (1992, 1995). Note that LRS (1993b) and RS (1992, 1994, 1995) define binary separation as 
the distance between the centers of mass of the binary components and, as mentioned above, we 
define it as the distance between the pressure (density) maxima of the components. For ease of 
comparison, the values of binary separation presented in this section that are based on our work 
do represent the separation between the centers of mass of our binary components. In this section 
and the rest of this manuscript, any binary separation which refers to the distance between the 
components centers of mass will be denoted as a cm . 

The analytical stability analyses of LRS (1993b), which utilized an approximate energy 
variational technique, also showed that simultaneous minima in E and J existed along sequences 
of constant mass, synchronized binaries with n = 0.5, 1.0 and 1.5. (According to LRS 1993b, the 
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minima occur at a cm = 2.99 on the n = 0.5 sequence and at a cm = 2.76 on the n = 1.0 sequence.) 
This analytical method cannot construct contact binaries and thus, according to our sequences, 
should not be able to identify minima on the n = 1.5 sequence. However, because the n at which 
the minima and points of contact coincide is ~ 2.0 in their study, LRS (1993b) do find minima for 
the n = 1.5 sequence at a cm = 2.55. 

In addition, approximate equilibrium sequences were constructed with the SPH techniques 
of RS (1992, 1995). The n = 1.0 sequence presented in LRS (1993b) has simultaneous minima 
at a cm = 2.90, which is closer to our value of a cm = 2.98 than the analytically determined 
a cm = 2.76 of LRS (1993b). The n = 1.5 sequence presented in RS (1995) has the minima at 
a cm = 2.67. No sequence with n = 0.5 has been published by RS. Table 2 contains a summary of 
the separations a cm of the models at the minima of the polytropic sequences as determined in LRS 
(1993b), RS (1995), and this work. For completeness, this table also gives the values of binary 
separation determined by this work in terms of the separation a between the pressure maxima of 
the components. 

Hachisu (1986b) also shows equilibrium sequences of n = 0.5 and 1.5 polytropes. However he 
presents his results as sequences of versus J 2 instead of E or J versus a. A comparison between 
Hachisu's (1986b) results and ours is given in Figure 3. To conform with Hachisu's notation, the 
quantities in this figure are normalized to 4ttG, Mt, and Vp, where Vt is the total volume of the 
binary system. Specifically, Qq has been divided by (AttGMt/Vt) and J 2 has been divided by 
(4ttGM t V^ /3 ). 

3.2. White Dwarf Sequences 

We have constructed equilibrium sequences for binary models with the zero-temperature WD 
equation of state. Because the HSCF technique requires that the maximum density of the desired 
model (which sets Mt) be given as input when this equation of state is used, it is impossible to 
build a single sequence that can be scaled to any desired Mt with this equation of state. Instead 
we have constructed separate WD sequences, each of which represents models with one specific 
value of Mt- We have constructed nine such sequences with Mt ranging from .298 Mq to 2.72 
Mq. Four representative sequences, with My = .500, 1.19, 2.03, and 2.72 Mq, are shown in Figure 
4. The other five sequences are displayed in Appendix A, along with the four WD sequences 
presented in Figure 4. (In addition, two WD sequences with Mt [near 2Mch] = 2.81 and 2.85 
are presented in Appendix A; they have been excluded from the discussion below because of their 
irregular nature.) The normalization in Figure 4 is the same as that in Figure 1. However, in 
this case Rs P h has been determined numerically by constructing a spherical white dwarf of mass 
M = Mt /2 in hydrostatic equilibrium. 

The separations a of the models at the points of contact, the minima, and the terminal 
points on the constructed white dwarf equilibrium sequences are shown in Figure 5 as a function 
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of the total binary system mass. Along the WD sequences, the separation of the terminal model 
gradually increases from a = 2.45 when M T = .298 M to a = 2.86 when M T = 2.72 M . As 
in the polytropic sequences presented in the previous section, simultaneous minima in E and 
J exist along each WD sequence. The point of contact on these sequences always occurs at a 
larger separation than the minima; the separation at which it occurs also gradually increases from 
a = 2.81 for M T = .298 M® to a = 3.05 for M T = 2.72 M Q (except for a slight decrease in this 
separation for the Mt = 1.19 M sequence). The separation of the model at the minimum of 
each sequence also increases from a = 2.70 for My = .298 M to a = 2.97 for Mt = 2.72 M ; 
however, in this case there is a somewhat more substantial decrease in this separation between the 
M T = .500 M & and M T = 2.36 M Q sequences (for which a = 2.73). 

3.2.1. Comparison with Previous Work 

Hachisu (1986b) has constructed double white dwarf binary sequences along which p m axj 
instead of Mt, was held constant. Figure 6 shows a comparison between the J versus Mt relations 
for the models with the minimum angular momentum on the Hachisu (1986b) sequences and those 
on our sequences. The angular momentum in this figure is normalized to 10 50 in cgs units and 
the mass is normalized to Mq. The comparison is excellent for low masses, however the relations 
deviate slightly for M T > 2 M . 

3.3. Neutron Star Sequences 

We have constructed equilibrium sequences for binary models with three realistic NS equations 
of state of varying compressibility (F, soft; FPS, medium; L, hard). As in the case of the WD 
equation of state, the desired Mt for each of these sequences must be specified prior to their 
construction. We have chosen to construct one sequence with Mt = 2.80 M for each of the three 
equations of state. These three sequences are displayed in Figure 7 and have been normalized to 
G, M, and Rs P h- The values of Rs p h were determined numerically by constructing a spherical 
star of mass M = Mr/2, in hydrostatic equilibrium, with each of the three NS equations of state. 

The terminal model occurs at the same separation on the F and FPS sequences (a = 2.59) 
and at a slightly wider separation (o = 2.62) on the L sequence. The locations of the minima in 
E are not very well defined on these sequences. For the F sequence, the E minimum likely occurs 
in the range 2.75 < a < 2.81; for the FPS sequence, in the range 2.72 < a < 2.81; and for the 
L sequence, in the range 2.72 < a < 2.84. The J minima are well defined and occur at nearly 
the same separation on all three sequences (on the F sequence at a = 2.80, on the FPS sequence 
at a = 2.78, and on the L sequence at a = 2.79). Like the white dwarf sequences, the minima 
always occur at a smaller separation than the points of contact. However, on all of the neutron 
star sequences these two points are very close together. The separation between points C and 
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M thus seems to be determined by the lower density regions of the objects since this is where 
the equations of state, which differ significantly only in the density regimes above nuclear density 
(2.8 x 10 14 gem -3 ), are similar. 

The scatter in E displayed near the minima of these sequences may result in part from our 
piecemeal reconstruction of Cook's (1995) tabular NS equations of state (see §2) and from the 
approximate manner in which we calculate the internal energy Ei n t of the models (we compute an 
effective polytropic index n e ff for each grid zone in the model and then use these spatially varying 
indices to calculate the internal energy according to Ei nt = J n e ffPdV). Because of these factors, 
we identify the J minima as the true minima along these sequences and, as before, have marked 
their location with the letter M. Another factor which must be kept in mind when studying the 
features of these NS sequences is that the virial error (VE), which measures the quality of the 
equilibria, is ~ 1(T 2 . This is one to two orders of magnitude higher than the VE of models along 
the polytropic and WD equation of state sequences and almost certainly results in part from the 
piecemeal forms of the NS equations of state used. Note that for the polytropic and WD equation 
of state, we typically very accurately pinpointed the location of the miminum by moving point B 
(one of the two points on the surface of the star that must be given as input to the HSCF code 
and whose position influences the separation of the binary components) one grid cell at a time in 
the region of the minimum in order to get as many models as possible with the grid resolution we 
used. However, since we had difficulty obtaining converged models for some portions of the NS 
sequences, we moved point B two grid cells at a time when constructing NS models. Thus the 
locations of the minima are slightly more uncertain for the NS equations of state than in the other 
studied equations of state. 

3.4. Nature of the Minima 

Turning points on equilibrium sequences, like the minima present in E and J along the 
sequences presented above, are usually associated with points of instability. The two types of 
instability which will be discussed below are secular instability and dynamical instability. For 
an instability to be classified as dynamical, according to the convention of LRS (1993b), it must 
conserve energy, angular momentum, and circulation. If an instability proceeds as a result of 
some mechanism that dissipates one of the conserved quantities, e.g., via viscosity or gravitational 
radiation, then it is classified as a secular instability (see also Tassoul [1978] for discussions of 
dynamical versus secular instabilities). A dynamical instability takes place on the dynamical 
timescale of the system; a secular instability takes place on the timescale of the particular 
dissipative mechanism which induces it. 

The minima in E and J on synchronous binary sequences have been identified as points of 
secular instability by previous authors. LRS (1993b) stated that the neighboring models adjacent 
to the model at the minimum of each sequence are uniformly rotating and therefore can only be 
reached with the aid of viscosity. Thus they concluded that this minimum cannot be associated 
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with a dynamical instability since viscosity does not preserve circulation. Their approximate 
analytical method predicts that on polytropic sequences with n < 0.7, dynamical instabilities exist 
at separations smaller than those of the minima. For the n = 0.5 sequence, these authors conclude 
that the dynamical instability sets in at a cm = 2.68. Hachisu (1986b and references therein) also 
labels turning points along his synchronous polytropic and WD sequences as points of secular 
instability. 

Because our hydrodynamics code does not include the dissipative effects of either the 
gravitational radiation reaction or the true fluid viscosity, we are unable to confirm the presence 
of a secular instability with a hydrodynamics simulation. However, it is possible to study the 
dynamical stability of a system with our code. In the following section, we present the results of 
our FDH tests of the dynamical stability of models on sequences selected from those presented 
above and compare these results with those of other authors. 

4. Hydrodynamic Tests of Stability 

To determine if a point of dynamical instability exists on an equilibrium sequence like those 
discussed in the previous section, the dynamical stability of individual models along the sequences 
may be tested with FDH simulations. We have done just that for the three polytropic sequences 
of §3.1; for a low mass (Mt = -500 M Q ) and a high mass (Mt = 2.03 M Q ) WD sequence from 
§3.2; and for two of the three realistic NS equations of state sequences of §3.3. 

All of these stability tests were performed in the rotating reference frame (§2.2) in order to 
minimize the dissipative effects of the numerical viscosity which results from the discrete nature 
of the computational simulation. The influence of numerical viscosity on the evolution of binary 
systems can be seen in Figure 8 which shows a comparison between simulations of one particular 
WD binary system (Mt = .500 M & and a = 2.63) performed in the inertial frame (dashed curve) 
and in the rotating frame (solid curve). This figure shows the evolution of the moment of inertia 
of the system /, as a function of time t. The evolution of / is more instructive than the evolution 
of the binary separation itself since the latter quantity can only be measured by discrete jumps in 
the location of the density maximum on the grid, whereas / varies smoothly with time. 

As can be seen from Figure 8, the binary appears to be dynamically unstable when the 
simulation is performed in the inertial frame (dashed curve), as / plummets on a timescale of 
2-3 Pi. By contrast, the same model evolved in the rotating frame is certainly not unstable to 
merger, as can be seen by the relatively steady behavior of its moment of inertia over time. Thus, 
simulations done in the rotating frame prevent the misidentification of models as dynamically 
unstable, which are so only because of numerical artifacts. In order to illustrate that the accuracy 
of the finite differencing scheme also influences the amount of numerical viscosity present in a 
simulation, Figure 8 also shows the same model evolved in the rotating frame but with a FDH code 
that used a first-order accurate advection scheme to compute the divergence terms in equations 



-14- 



(1-4) (dot-dashed curve). This figure indicates that the accuracy of the code has an effect on the 
evolution of this system that is similar to that of the flow of the fluid through the grid in the 
inertial frame. 

Models along the equilibrium sequences presented in §3 were constructed with a grid resolution 
of 128 x 128 x 128. However, a FDH simulation with this resolution cannot be done on the MP-1 
computer at LSU. So portions of the selected sequences mentioned above have been recomputed 
on 64 x 64 x 64 grids. It is the stability of models on these new sequences which actually has 
been tested. For completeness, Table 3 gives the separations of the points of contact, minima, 
and terminal models for both the 64 3 and the 128 3 versions of the polytropic and WD sequences 
discussed below. (We do not show this comparison for the NS sequences because we have not 
determined or were unable to accurately determine the location of some of these points on the 64 3 
NS sequences.) 

4.1. White Dwarf Sequences 

4.1.1. M T = 2.03 M Q 

In the lower panel of Figure 9, I(t) is shown for WD binaries (Mt = 2.03 M & ) with separations 
ranging from a = 2.80 (triple dot-dashed curve, a dumbbell model just past the point of contact) 
to a = 2.53 (solid curve; the terminal model on the sequence). For the sake of convenience, the 
relevant equilibrium sequence itself is reprinted in the upper panel of Figure 9. The moments 
of inertia of the binary models at the points of contact, minima, and termination along the 
equilibrium sequence (constructed on a 128 x 128 x 128 grid) are labelled with the letters M, C, 
and T, respectively. Note that because the binary models used in the FDH stability tests were 
from sequences constructed on 64 x 64 x 64 grids, there may be a slight offset between the initial 
values of / for these models and the points marked as M, C, and T on the I(t) plot since they 
correspond to the 128 x 128 x 128 sequence. 

All of the models tested on the sequence are dynamically stable. Thus we conclude that 
no point of dynamical instability exists along this sequence. In Figure 10, we present isodensity 
images of an example of the evolution of a stable binary. The images are from the simulation of 
the model at the minimum of the Mt = 2.03 M sequence, corresponding to the dot-dashed curve 
of Figure 9. Since the simulation was performed in the rotating frame, the binary does not appear 
to pivot very much in this set of images even though the simulation was carried out for ~ \Pi . 

4.1.2. M T = . 500 M 

Our tests of the Mt = 2.03 Mq sequence began with models of wider separation and continued 
up the sequence to the terminal model. Because no point of dynamical instability was found on 
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this sequence, we performed our tests of the Mr = .500 M & sequence on models much closer to 
the terminal model, knowing that we could work our way back down the sequence to the point of 
dynamical instability if we found models which were unstable. (This methodology assumes that 
all models which have separations less than that of the model located at the point of dynamical 
instability will also be unstable.) But, as can be seen from the bottom panel of Figure 11, both a 
model at a = 2.63 (dashed curve) and the terminal model at a = 2.45 (solid curve) are dynamically 
stable against merger. Thus we again infer that there is no point of dynamical instability along 
this sequence. 

4.2. Polytropic Sequences 

4.2.1. n = 1.5 

Based on the results of the stability tests of the WD sequences, we began our investigation 
of this sequence with the terminal model (a = 2.45). Because an n = 1.5 polytrope is supposed 
to be a fair representation of a low mass WD, we have plotted the I(t) of the terminal model 
(a = 2.45) in Figure 11 (dot-dashed curve) along with the models of the low mass WD sequence. 
The polytropic model has been scaled to represent a binary with Mt = .500 M & and a Rsph equal 
to that of a spherical WD with the same mass. This sets K according to equation 13. As can be 
seen, the evolution of the n = 1.5 terminal model closely follows that of the low mass WD terminal 
model. From its stability, we infer that the remainder of the models on the n = 1.5 sequence are 
also stable and thus no point of dynamical instability exists along this sequence either. This result 
conflicts with the SPH simulations of RS (1995) which identified a point of dynamical instability 
at a cm ~ 2.4 on their n = 1.5 sequence. 

4.2.2. n = 1.0 

As Figure 12 illustrates, we have performed dynamical stability tests of four binaries on the 
n = 1.0 sequence with separations ranging from a = 3.28 (triple dot-dashed curve) to a = 2.62 
(solid curve; terminal model). All were stable against merger on a dynamical timescale. The 
moment of inertia of the terminal model actually increases by about 20% over the timescale 
t = APj depicted here; the moment of inertia of the initially most widely separated model 
(a = 3.28) decreases by about the same percentage. The two other models exhibit behavior in 
between these two extremes. Note that after an initial dip, the moment of inertia of the model at 
the minimum of the sequence (o = 2.88, the dot-dashed curve) is nearly constant. These results 
conflict with the SPH simulations of RS (1992) which identified a model with being 
dynamically unstable. 
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4-2.8. n = 0.5 

As Figure 13 illustrates, we have tested the dynamical stability of five binaries on the n = 0.5 
sequence with separations ranging from a = 3.41 (long dashed curve) to a = 2.77 (solid curve; 
terminal model). Both the model at the minimum of the sequence (a = 3.10-short dashed curve) 
and the terminal model were unstable to merger on a dynamical timescale. The other three more 
widely separated models, including the model (a = 3.17-dot-dashed curve) located just prior to 
the minimum of the sequence, were stable. Thus dynamical instability sets in at the minimum 
of this sequence (a = 3.10, a cm = 3.20). The SPH simulations of RS (1994) identified a point of 
dynamical instability along this sequence at a cm = 2.97. 

4.3. Neutron Star Sequences 

4-3.1. Equation of State L 

As Figure 14 illustrates, we have tested the dynamical stability of a NS binary system near the 
minimum (a = 2.79) and a more widely separated model (a = 3.26) on the realistic neutron star 
equation of state L sequence. The model near the minimum of the sequence became unstable to 
dynamical merger in one Pi. The more widely separated model appears to be stable; its behavior 
resembles that of the most widely separated models tested on the n = 1.0 and 0.5 sequences 
(see Figures 12 and 13). Although we have not accurately identified the separation at which the 
dynamical instability sets in along this sequence, it has obviously already done so at a = 2.79. By 
analogy with the n = 0.5 sequence, it is possible to infer that the onset of dynamical instability is 
associated with the region of the minimum energy and angular momentum along the sequence. 

4-3.2. Equation of State F 

As Figure 15 illustrates, we have tested the dynamical stability of the model at the minimum 
(a = 2.80) of the realistic neutron star equation of state F sequence. Like the model near the 
minimum of the equation of state L sequence, the model at the minimum of this sequence is also 
unstable to dynamical merger. Although we have not tested the stability of a widely separated 
model on this sequence, we infer, by analogy with the n = 0.5 and equation of state L sequences, 
that there exist more widely separated models on this sequence that are stable and that the point 
of onset of the dynamical instability may be associated with the region of the minimum. 
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4-3.3. Equation of State FPS 

We have not tested the stability of any of the models on the realistic NS equation of state 
FPS sequence. However, since our tests of the stiff equation of state L and the soft equation of 
state F produced similar results, we expect that the stability properties of the models on this 
equation of state of medium stiffness will resemble the properties of models on these other realistic 
NS equations of state sequences. 

5. Discussion and Conclusions 

We have examined the dynamical stability of synchronously rotating, equal-mass binaries 
with polytropic, zero-temperature white dwarf, and realistic neutron star equations of state. 
Specifically, we tested the dynamical stability of individual models constructed along equilibrium 
sequences of binaries with the same total mass My and equation of state but decreasing separation, 
in order to determine if any models on these sequences were unstable to merger on a dynamical 
timescale. 

Our stability analyses started with two white dwarf (WD) equation of state sequences, one 
with a low total mass (Mt = -500 M ) and one with a fairly high total mass (Mr = 2.03 M ), 
used as representatives of the nine regular WD equation of state sequences we constructed (which 
ranged in mass from My = .298 to 2.72 Mq). Our simulations indicate that no points of dynamical 
instability exist on either of these two sequences. We have inferred from this result, that it is 
likely (although not certain) that the other WD sequences are also dynamically stable. This being 
the case, we expect that WD mergers will only happen via secular processes and that it will not 
be possible to properly simulate the coalescence of equal-mass double white dwarf binaries using 
explicit FDH (or SPH) techniques. It still remains to be seen whether or not double white dwarf 
binaries having unequal mass components are susceptible to merger on a dynamical timescale. 

Our examination of the n = 1.5 and n = 1.0 polytropic sequences also identified no points 
of dynamical instability. This result conflicts with the published results of RS (1992, 1995) who 
identified a dynamically unstable binary between the minimum and the terminal point of the 
n = 1.5 sequence and another just past the minimum of the n = 1.0 sequence. However, like RS 
(1994), our test of the n = 0.5 sequence does indicate the presence of a dynamical instability. In 
our simulations, this instability sets in at the minimum of the sequence. RS (1994) locate the point 
of dynamical instability at a cm = 2.97. Although they do not state where the minimum of their 
sequence is located, the analytical work of LRS (1993b) places it at ct crn = 2.99. Recall that LRS 
(1993b) label this minimum as a secular instability and they predict that a dynamical instability 
sets in at a smaller separation (a cm = 2.68) on this sequence. However, in light of the simulations 
of RS (1994) and of this work, it seems likely that the minimum itself may be associated with the 
onset of the dynamical instability. 
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The differences between our results and those of RS (1992, 1995) for the sequences with the 
softer equations of state may result from the fact that SPH has difficulty modeling low density 
regions that are more extensive in systems with softer equations of state. It is interesting to note 
that the analytical work of LRS (1993b) predicts that points of dynamical instability exist only 
on sequences for which n is less than « 0.7. It would be intriguing to test the stability of other 
polytropic sequences with our FDH code in order to determine at what value of n the dynamical 
instability first appears. Another point that is worth emphasizing again, is that our stability tests 
were run in the rotating frame of reference. We believe that we would have misidentified some 
truly stable models as being unstable had these tests been performed in the inertial frame (see the 
discussion in §4 associated with Figure 8). 

Our stability analyses of the realistic neutron star (NS) equations of state sequences identified 
models at (or near) the minima of the soft equation of state F and the hard equation of state 
L as being dynamically unstable. However, because of computing resource constraints we have 
been unable to determine whether or not the onset of the dynamical instability takes place in the 
regions of the minima of these sequences. If the n = 0.5 sequence is taken as an example, one 
might infer this to be the case. Although we did not test the stability of the medium equation of 
state FPS sequence we constructed, we expect that such tests would yield results similiar to those 
of the other two NS equations of state. If further simulations of the n = 0.5 and NS equations 
of state sequences do identify the minimum as the point of dynamical instability, the question 
will arise as to why this turning point does not also mark the onset of dynamical instability on 
the sequences with softer equations of state. At this time we are unable to provide a physical 
explanation for this possibility (see, however, the discussion in LRS [1993b]). 

All of the binary models that we have identified as stable against dynamical merger in §4 
have not necessarily become steady-state configurations by the end of our simulations. In fact, at 
the end of our simulations the moments of inertia of some of the stable models are still decreasing 
gradually and those of others are still increasing gradually. However, these binaries did not become 
unstable to merger on a dynamical timescale (i.e., ~ a few Pj). 

It is interesting to note that this gradual (secular-type) evolution of our models is such that 
models before the minimum of a sequence (i.e., models with separations larger than that of the 
model at the minimum) tended to exhibit decreasing moments of inertia and those after the 
minimum tended to exhibit increasing moments of inertia. (This trend appears in the simulations 
of models on the low mass white dwarf sequence and is more pronounced in the tests of models 
with stiffer equations of state.) LRS (1994b) have indicated that viscosity does indeed cause the 
orbits of models before the minimum of a sequence to decay. They also suggest that as a result 
of the action of viscosity, the orbit of a binary past the minimum ". . . can either expand ... as the 
system is driven to a lower energy, stable synchronized state, or it can decay ... as the stars are 
driven to coalescence." 

The LRS (1994b) description of the action of viscosity is consistent with our results if the 
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gradual evolution of the moments of inertia of our models is attributed to the effects of the 
remaining numerical viscosity in our code (which acted in our case to increase the moments of 
inertia of dynamically stable models past the minima). By carrying out simulations in the rotating 
reference frame (as is illustrated in Figure 8), we have overcome the primary effect of numerical 
viscosity, which is to dissipate orbital motions and drive systems to coalescence. In simulations 
performed in the rotating frame, it may be a higher order effect of numerical viscosity on any 
internal motions that develop in the systems that causes the observed secular-type evolution of 
our binary models. Note that this higher order effect of numerical viscosity cannot be identified 
as the mechanism responsible for driving any of our unstable models to coalescence because 
the timescale of their evolution was dynamical and not gradual. In addition, had this aspect of 
numerical viscosity played a significant role in the evolution of the terminal model on the n = 0.5 
sequence, its moment of inertia should have increased gradually, but instead its moment of inertia 
decreased and the model merged on a dynamical timescale. 

Even though our work indicates that there is no point of dynamical instability along the 
n = 1.0 sequences, it would still be possible for an explicit hydrodynamics code to follow the 
merger of a close binary with this equation of state if it was assumed to represent a neutron 
star-neutron star system. This is because the timescale, r, on which gravitational radiation 
would drive a close neutron star-neutron star binary to coalescence is on the order of its initial 
orbital period, Pj. According to Shapiro & Teukolsky (1983), a point mass approximation to this 
timescale is 

T 256G 3 2M 3 ' [ ' 

If the binary at the minimum of our n = 1.0 sequence is assumed to have components with 
M = 1.4 M and Rs P h = 10 km, then r = 2.5 ms = 1.6 Pj. (Note that for the binary at the 
minimum of our Mt = 2.03 M WD equation of state sequence, rPf 1 = 2.2 x 10 7 and for the 
binary with a = 2.63 on our My = .500 M & WD equation of state sequence, rPf 1 = 6.4 x 10 9 .) 
Hence, if the effects of the gravitational radiation reaction were accounted for, as other authors 
have done (e.g., Oohara & Nakamura 1990; Davies et al. 1994; Zhuge, Centrella, & McMillan 1994; 
Ruffert, Janka, & Schafer 1996), the merger of the n = 1.0 binary would proceed on a timescale 
comparable to the dynamical timescale. 

The variation between the results of our stability analyses of binaries with the softer equations 
of state and those of RS (1992, 1995) emphasizes the importance of comparing the results of 
hydrodynamic simulations performed with different numerical techniques. Certainly the inclusion 
of more complex physics, such as the effects of the gravitational radiation reaction and of full 
general relativity should be the aim of research in this field. However, an effort should be made to 
reach agreement on the results of much simpler, Newtonian simulations in order that their results 
and the results of more complex simulations can be confidently presented as guides to those who 
will be designing and building gravitational radiation detectors and interpreting the data collected 
at future gravitational radiation observatories. 
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A. Equilibrium Sequences 



Eleven WD equation of state equilibrium binary sequences, each with a different My, are 
displayed in Figure 16. See §3.2 for details. 
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FIGURE CAPTIONS 

Fig. l.-Polytropic Equilibrium Sequences. Sequences of binaries with polytropic indices 
n = 0.5, 1.0, and 1.5 are displayed. Each sequence displays the total energy E and the total 
angular momentum J of synchronized equilibrium binaries with the same total mass Mt but 
changing separation a. See the text (§3.1) for details on how to scale the Mt of each polytropic 
sequence. Quantities have been normalized to G, M = Mr/2, and Rs P h, where Rs P h is the radius 
of a spherical polytrope of mass M and polytropic index n. Points M mark binary systems with 
the minimum E and J along each sequence. Points C mark the system where the surfaces of the 
stars come into contact. Points T mark the terminal model. 

Fig. 2.-Images of a Detached Binary and a Dumbbell. Example isodensity images of a 
detached binary and a contact binary, or "dumbbell." These two binaries have n = 1.0; the 
separation of the detached binary shown in Figure 2a is a = 3.28; the separation of the dumbell 
shown in Figure 2b is a = 2.70. The density level is 5.0 x 10~ 3 of the maximum density. 

Fig. 3.-Comparison with Hachisu's Polytropic Sequences. The square of the angular velocity 
fio versus the square of the total angular momentum J of equilibrium binaries with the same 
total mass Mt but decreasing separation a is shown for both the n = 0.5 and 1.5 sequences of 
Hachisu (1986b, individual models connected by solid lines) and for our sequences (individual 
models marked with +'s). The quantities are normalized to AirG, Mt, and Vr, the total volume 
of the system. 

Fig. 4.-White Dwarf Equilibrium Sequences. The same as Figure 1 but for binaries with 
the zero-temperature white dwarf equation of state. Because, unlike the polytropes, the Mt of 
these systems cannot be scaled, separate sequences must be constructed for binaries with different 
Mt- Four representative sequences with Mt = .500, 1.19, 2.03, and 2.72 M & are shown here. In 
Appendix A, sequences for seven other values of Mt are shown, along with those given here. The 
normalization is the same as in Figure 1, but here Rs P h is the radius of a spherical model with the 
WD equation of state and mass M = My/2. 

Fig. 5.-Separation of Models on White Dwarf Equilibrium Sequences. The separation a of 
the models at the points of contact (x), the minima (asterisks), and the terminal points (plusses) 
on the white dwarf equilibrium sequences are shown as a function of the total mass Mt of the 
binaries on those sequences. 

Fig. 6.-Comparison with Hachisu's White Dwarf Sequences. The total angular momentum 
J versus the total mass Mt of the models with the minimum J on each of Hachisu's (1986b) 
constant maximum density sequences (individual models marked with x's) and on each of our 
constant Mt sequences (individual models marked with +'s). The angular momentum has been 
normalized to 10 50 in cgs units and the total mass to Mq. 

Fig. 7. -Neutron Star Equilibrium Sequences. The same as Figure 1 but for binaries with the 
F, FPS, and L realistic neutron star equations of state. The My on each of these sequences is 2.80 
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M (unlike polytropes, the Mt of these systems cannot be scaled). The normalization is the same 
as in Figure 1 but here Rs P h is the radius of a spherical model constructed with the particular NS 
equation of state of the sequence and a mass M = My/2. 

Fig. 8. -Comparison of Stability Tests Performed with Different Numerical Techniques. The 
moment of inertia /, normalized to (MPJ| . ), as a function of time t, normalized to the initial 
orbital period Pj, is given for dynamical stability tests of a binary with the zero-temperature 
white dwarf equation of state, Mt = 0.500 M Q , and a = 2.63. The solid curve shows the results 
of a simulation performed with our second-order accurate FDH code in the rotating frame; the 
dashed curve shows the test performed with the same code but in the inertial frame; and the 
dot-dashed curve shows the test performed in the rotating frame, but with the advection terms in 
the FDH code reduced to first-order accuracy. 

Fig. 9.-Stability Tests of the M T = 2.03 M White Dwarf Sequence. The variables in 
the lower plot are the same as Figure 8, but the simulations all are performed in the rotating 
reference frame with the second-order accurate FDH code. The solid curve is the stability test 
of the terminal model on the sequence (a = 2.53); the dashed curve is the test of the model 
with a = 2.63; the dot-dashed curve is the test of the model at the minimum (a = 2.67); and 
the triple dot-dashed curve is the test of the model just past point of contact (a = 2.80). The 
moments of inertia of the binary models at the points of contact, minima, and termination along 
the equilibrium sequence are labelled with the letters C, M, and T, respectively. For convenience, 
the sequence itself has been reprinted from Figure 4 in the upper plot. 

Fig. 10. -Evolution of Stable Binary. Isodensity images of a stable binary with the WD 
equation of state, Mt = 2.03 Mq, and a = 2.67. The density level is 5 x 10 -4 of the maximum 
density. Figure 10a was taken at t = 0.0 Pr, 10b at t = 0.8 Pj; 10c at t = 1.5 Pr, 10d at t = 2.3 Pr, 
lOe at t = 3.0 Pr; and lOf at t = 3.8 Pj. Since the simulation was performed in the rotating frame, 
the binary does not appear to pivot very much in this set of images. 

Fig. ll.-Stability Tests of the M T = .500 M White Dwarf and the n = 1.5 Sequences. The 
same variables and code as used in Figure 9. The solid curve is the simulation of the terminal 
model (a = 2.45) on the Mt = -500 Mq white dwarf sequence (the equilibrium sequence itself is 
reprinted from Figure 4 in the upper plot); the dashed curve is the test of the model on this same 
sequence with a = 2.63; and the dot-dashed curve is the test of the terminal model (a = 2.45) on 
the n = 1.5 polytropic sequence. 

Fig. 12. -Stability Tests of n = 1.0 Sequence. The same variables and code as described in 
Figure 9. The solid curve is the test of the terminal model (a = 2.62); the dashed curve is the test 
of the model with a = 2.70; the dot-dashed curve is the test of the model at the minimum of the 
sequence (a = 2.88); and the triple dot-dashed curve is the test of the model with a = 3.28. The 
equilibrium sequence itself is reprinted from Figure 1 in the upper plot. 

Fig. 13.-Stability Tests of n = 0.5 Sequence. The same variables and code as described in 
Figure 9. The solid curve is the test of the terminal model (a = 2.77); the short dashed curve is 
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the test of the model at the minimum of the sequence (a = 3.10); the dot-dashed curve is the test 
of the model with a = 3.17; the triple dot-dashed curve is the test of the model with a = 3.24; 
and the long dashed curve is the test of the model with a = 3.41. The equilibrium sequence itself 
is reprinted from Figure 1 in the upper plot. 

Fig. 14. -Stability Tests of Equation of State L Sequence. The same variables and code as 
described in Figure 9. The solid curve is the test of the model near the minimum of the sequence 
(a = 2.79); the dashed curve is the test of the model with a = 3.26. The equilibrium sequence 
itself is reprinted from Figure 7 in the upper plot. 

Fig. 15.-Stability Tests of Equation of State F Sequence. The same variables and code as 
described in Figure 9. The solid curve is the test of the system at the minimum of the sequence 
(a = 2.80). The equilibrium sequence itself is reprinted from Figure 7 in the upper plot. 

Fig. 16.-White Dwarf Equilibrium Sequences. M T = .298, .500, .803, 1.19, 1.63, 2.03, 2.36, 
2.58, 2.72, 2.81, and 2.85 M . 



Table 1: Lane-Emden Constants, m n and 



n m n r n 

0.5 3.7871 2.7528 

1.0 3.14159 3.14159 

1.5 2.71406 3.65375 



"From Chandrasekhar (1967) 



Table 2: Separation of Models at Minima of Polytropic Sequences 



Technique, Authors 
Analytic, LRS 
SPH, RS 
HSCF, New & Tohline 
HSCF, New & Tohline 



n = 0.5 n = 1.0 n = 1.5 

2.99 a 2.76 a 2.55 a 

2.90 a 2.67 b 

3.20 c 2.98 c 2.77 c 

3.11 d 2.89 d 2.70 d 



"From LRS (1993b), separation between components' centers of mass 
''From RS (1995), separation between components' centers of mass 
"From this work, separation between components' centers of mass 
d From this work, separation between components' pressure maxima 
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Table 3: Comparison of 64 3 and 128 3 Sequences 



EOS, M T 
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a(T) 


WD, 2.03 M 
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2.94 
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2.94 


3.11 


2.76 
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